Comprehensive prediction of immune microenvironment and hot and cold tumor differentiation in cutaneous melanoma based on necroptosis-related lncRNA

As per research, causing cancer cells to necroptosis might be used as a therapy to combat cancer drug susceptibility. Long non-coding RNA (lncRNA) modulates the necroptosis process in Skin Cutaneous Melanoma (SKCM), even though the precise mechanism by which it does so has yet been unknown. RNA sequencing and clinical evidence of SKCM patients were accessed from The Cancer Genome Atlas database, and normal skin tissue sequencing data was available from the Genotype-Tissue Expression database. Person correlation analysis, differential screening, and univariate Cox regression were successively utilized to identify necroptosis-related hub lncRNAs. Following this, we adopt the least absolute shrinkage and selection operator regression analysis to construct a risk model. The model was evaluated on various clinical characteristics using many integrated approaches to ensure it generated accurate predictions. Through risk score comparisons and consistent cluster analysis, SKCM patients were sorted either high-risk or low-risk subgroups as well as distinct clusters. Finally, the effect of immune microenvironment, m7G methylation, and viable anti-cancer drugs in risk groups and potential clusters was evaluated in further detail. Included USP30-AS1, LINC01711, LINC00520, NRIR, BASP1-AS1, and LINC02178, the 6 necroptosis-related hub lncRNAs were utilized to construct a novel prediction model with excellent accuracy and sensitivity, which was not influenced by confounding clinical factors. Immune-related, necroptosis, and apoptosis pathways were enhanced in the model structure, as shown by Gene Set Enrichment Analysis findings. TME score, immune factors, immune checkpoint-related genes, m7G methylation-related genes, and anti-cancer drug sensitivity differed significantly between the high-risk and low-risk groups. Cluster 2 was identified as a hot tumor with a better immune response and therapeutic effect. Our study may provide potential biomarkers for predicting prognosis in SKCM and provide personalized clinical therapy for patients based on hot and cold tumor classification.

Necroptosis-related lncRNAs were detected by Pearson correlation analysis based on coefficients > 0.4 and P < 0.001. Next, differential expression necroptosis-related lncRNAs were identified using the limma R package, in which the screening standard is |log2 fold change|> 2 and false discovery rate (FDR) < 0.05.
Construction of the risk model. Merge clinical information from tumor samples with all differentially expressed lncRNAs data. Overall prognostic impact of necroptosis-associated lncRNAs was evaluated using univariate Cox (uni-Cox) regression analysis to filter potential lncRNAs (p < 0.05). Next, a least absolute shrinkage and selection operator (LASSO) regression analysis was run on candidate necroptosis-associated lncRNAs to identify hub lncRNAs and to create a prediction model through multivariate Cox regression. 1000 cycles of LASSO regression analysis were performed using tenfold cross-validation at a significance level of 0.05. Each cycle was randomly grouped 50 times, dividing the sample equally into two groups each time. Furthermore, the LASSO regression results were output only when the error was minimum. The risk score formula is as described below: where exp(lncRNA) is the relative expression of signature lncRNAs, and coef (lncRNA) is the weighted correlation coefficient of signature lncRNAs. Risk classification of SKCM patients by comparing the risk score of the patient with the median score.
Constructing and evaluating the nomogram. Univariate cox (uni-Cox) and multiple cox (multi-Cox) regression analyses were performed to assess the relationship between risk score and clinical features. And based on these results, a nomogram was created to predict patient outcomes using the RMS R package. Furthermore, the correction curves were constructed based on the Hosmer-Lemeshow test to verify the predictive reliability regarding the nomogram. The model's 1-, 3-and 5-year time-dependent receiver operating characteristics (ROC) curves were built by the timeROC R package. Finally, we used decision curve analysis (DCA) to further assess the predictive accuracy of this nomogram. www.nature.com/scientificreports/ Gene set enrichment analysis. Determine if there are significant differences throughout the functional pathways between high-and low-risk subgroups, the GSEA 4.2 software was used to discover the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways notably enriched in the two risk groupings (the criterion was defined as FDR < 0.05).
Correlation analysis of immunologic characteristics and m7G methylation. Immune infiltration by different algorithms was downloaded from TIMER2.0 (http:// timer. cistr ome. org/) and merged with the risk score. Then, the association between the risk model and immune infiltration was estimated derived from Spearman's rank correlation coefficient. The two risk subgroups were also compared using a single-sample gene set enrichment analysis (ssGSEA) to assess whether or not there was a massive distinction in immune function and immune cell enrichment. Moreover, the high-and low-risk group were compared in relation to the tumor microenvironment (TME) score, and expression of immune checkpoint-related genes. To explore the relationship between necroptosis-related lncRNAs and m7G methylation, we compared the expression of m7g-related genes between high-and low-risk group segmented according to this risk signature.
Estimation of anticancer drug sensitivity for SKCM. The "pRRophetic" R package, founded by Genomics of Drug Sensitivity in Cancer (GDSC)(https:// www. cance rrxge ne. org/), was employed to forecast the difference in drug sensitivity for both risk groups.
Consensus cluster analysis based on the hub lncRNAs. The ConsensuClusterPlus R package was implemented to separate SKCM patients into potential subgroups depending on the expression of hub lncRNAs by the k-means clustering algorithm. We set the maximum k value as 9 and determine the optimal K value through the combination of CDF and consensus matrices. Through using Rtsne and survival R package, we examined the subtype distribution and assessed the classification effectiveness including using t-distributed stochastic neighbor embedding (t-SNE), principal component analysis (PCA), and Kaplan-Meier survival analysis. Moreover, the immune status, m7G methylation-related genes, and clinical potential chemotherapeutic drug clusters were also analyzed simultaneously.

Results
Identify necroptosis-related lncRNAs in SKCM. There are 386 necroptosis-related lncRNAs identified from the data of TCGA and GTEx, as the standard is the coefficients > 0.4 and P < 0.001. After that, flowing the differential expression analysis, 87 necroptosis-related lncRNAs were found to display significantly differential expression with the screen value as |logFC |> 2 and FDR < 0.05. Among the expression profile of these lncRNAs in tumor samples was 57 upregulated and 30 downregulated. The network map, heat map, and volcano map are shown in Fig. 1.

Construction and verification of prognosis risk model.
On the basis of uni-cox regression analysis, 33 potential necroptosis-related lncRNAs that are substantially linked with overall survival were reported (p < 0.05) ( Fig. 2A). we also produced a heat map which shows how these lncRNAs were conveyed in both cancer and normal samples (Fig. 2B). Then, to construct the risk model, 6 hub necroptosis-associated lncR-NAs were selected through LASSO regression analysis. Furthermore, we visualized the results of Lasso regression and cross-validation (Fig. 2C,D). Among these hub lncRNAs, only the relationship between LINC01711 with BCL2 and AXL was negative, while the rest were positive (Fig. 2E). The risk score was calculated by formula: USP30-AS1 × (− 0.6073) + LINC01711 × (− 0.3162) + LINC00520 × (0.2742) + NRIR × (− 0.6839) + BASP1-AS1 × (0.3524) + LINC02178 × (1.0015). The results of scatter plots and risk curves, which were consistent across the training, testing, and overall set, show that patients have shorter survival times as risk scores increase ( Fig. 3A-F). USP30-AS1, LINC01711, and NRIR were clearly downregulated in the high-risk group, while LINC00520, BASP1-AS1, and LINC02178 presented opposite results in the heatmaps ( Fig. 3G-I). The K-M curves were consistent across the training, test, and total pools, showing that patients with low-risk scores had significantly better overall survival (OS) than the high-risk group, implying a worse prognosis for high-risk patients ( Fig. 3J-L). Furthermore, the K-M curves' results among the tumor's different stages were consistent with the front results ( Fig. 3M,N).

Construction of nomogram.
Compared with age, gender, and cancer stage, the hazard ratio (HR) of the risk score in uni-Cox and multi-Cox regression were 1.269 and 1.257 (p < 0.05) (Fig. 4A,B). Therefore, this risk signature may act as stable and specific biomarkers to predict the OS of melanoma patients. Upon the basis of risk score, age, and cancer stage, a nomogram was constructed to estimating the 1-, 3-and 5-year OS occurrences of SKCM patients (all p < 0.05 in multi-Cox) (Fig. 4C). Furthermore, the calibration curve result shows that the nomogram with an excellent concordance for 1-, 3-and 5-year (Fig. 4D). Additionally, by comparing the ROC curves of risk signature with age, gender, and cancer grade at 1-, 3-and 5-year, where the risk scores at 1-,3-and 5-year with AUC values of 0.704, 0.685 and 0.715 respectively, were significantly better than the other three clinical characteristics ( Fig. 5A-C). We also constructed DCA curves for the nomogram at 1-, 3-, 5-year compared to age and cancer grade, all of which showed that the risk profile had a more accurate predictive power ( Fig. 5D-F). These findings suggest that the Nomogram derived from this risk signature can accurately predict the prognosis of melanoma patients over the next 5 years. www.nature.com/scientificreports/ GSEA. The KEGG pathway enrichment analysis was carried out in the overall set to analyze the biological function of the risk signature. The top 10 pathways increased in the low-risk group were mainly associated with immunity function, including those involved in B cell, T cell, Toll-like receptor signaling, and natural killer cell-mediated cytotoxicity. Besides, the rest of the pathways were concerned with apoptosis, Cytokine-cytokine receptor interaction, and inflammation-induced immune response. Incredibly, there was no significant pathways enrichment in high-risk groups (FDR < 0.05) (Fig. 6).

Investigation of immune factors and m7G methylation-associated genes in risk models. From
the results on the bubble plot, the vast majority of immune cells were correlated with the low-risk group. The high-risk group was associated with a few immune cells, such as T cell NK and Mast cell at XCELL, Myeloid dendritic cell at QUANTISEQ, and Macrophage M1, M2 at CIBERSORT (Fig. 7A). For CIBERSORT platform, immune cell infiltration prediction showed that T cell CD8+, T cell CD4+ memory activated, Macrophage M1, NK cell activated, T cell gamma delta, and B cell memory was negative with risk scores. B cell naive, Myeloid dendritic cell activated, T cell CD4+ naive, Mast cell activated, Eosinophil, Macrophage M0, Macrophage M2, and NK cell resting were positive with risk scores (P < 0.05) (Fig. 8). The TME scores (containing stromal, immune, and ESTIMATE scores) were substantially higher in the low-risk group (P < 0.05) (Fig. 7B-D). We further observe that compared to the high-risk group, the low-risk group demonstrated significantly elevated levels of 15 immune cells score and 12 immune functions (***P < 0.001, **P < 0.01) (Fig. 7E,F). Only the scores of mast cell and type II IFN responses did not differ in the two groups. Above of this suggests that the low-risk www.nature.com/scientificreports/ group appears to be more strongly associated with the immune microenvironment. Furthermore, by differential expression analysis, the two groups showed significant differences in the expression profiles of 44 immunerelated checkpoint genes (***P < 0.001, *P < 0.05) (Fig. 7G). Among them, TNFRSF14, CD276, VTCN1, and CD44 were enriched in the high-risk group, whereas the rest genes were concentrated in the low-risk group. Therefore, we can use risk grouping to rationalize the selection of immune checkpoint inhibitors. Between risk groups, there were statistically significant changes in the expression levels of 17 genes involved in m7G methylation. In comparison to the low-risk subgroup, the high-risk subgroup exhibited markedly lower amounts of DCP2, IFIT5, NCBP2, and EIF4E3, while higher levels of EIF3D, EIF4A1, NSUN2, GEMIN5, AGO2, EIF4E, NCBP3, WDR4, LARP1, SNUPN, NUDT3, DCPS, and METTL1 (***P < 0.001, **P < 0.01, *P < 0.05) (Fig. 7H). www.nature.com/scientificreports/ Relevance between risk model and drug sensitivity. Necroptosis-related hub lncRNAs correlate with several common melanoma therapeutic agents. Low-risk patients were more responsive to Axitinib, Bosutinib, Cisplatin, Dasatinib, Gefitinib, Nilotinib, and Sunitinib than than high-risk patients. In contrast, the finding for imatinib was polar opposite (P < 0.05) (Fig. 9). Therefore, individualized selection of sensitive drugs for the risk subgroup was expected to improve the effectiveness of clinical treatment.

Differentiate between cold and hot tumors.
According to the results shown in the CDF diagram, consensus matrix plots, and fragment plots, the optimal K value was 2 (Fig. 10A,B). Then, we reclassified SKCM patients into two clusters depending on the expression of hub lncRNAs. The results of PCA and t-SNE show that the distribution of samples in clusters does not exactly match the risk subgroups (Fig. 10C-F). The Sankey plots show the relationship between the two grouping methods, where Cluster1 strongly linked to the highrisk group and Cluster 2 to the low-risk group (Fig. 10G). Following the K-M analysis, Cluster 2 exhibited a superior survival outcome (P = 0.002) (Fig. 10H). By evaluating the immune microenvironment for clusters, we found that Cluster 2 had a higher degree of immune infiltration and a significantly higher TME score than Cluster 1 (P < 0.05) (Fig. 11A-D). Only TNFRSF14, CD276, and VTCN1 were significantly expressed in Cluster 1. Conversely, most immune checkpoint genes were significantly activated in Cluster 2 (***P < 0.001; **P < 0.01; *P < 0.05) (Fig. 11E), especially the HAVCR2, LAG3, and CD274 genes, which were thought to contribute to important functions in hot tumors. Accordingly, cluster 2 was considered the hot tumor. The expression levels of m7G methylation-related genes DCP2, IFIT5, AGO2, NCBP2, EIF4G3, NUDT4, and EIF4E3 were higher, while NSUN2, EIF4E, WDR4, DCPS, and METTL1 were lower in cluster 2 in comparison with cluster 1 (Fig. 11F).
We also analyzed the differences in drug IC50 between clusters and showed that AKT.inhibitor.VIII, Axitinib, Cisplatin, Dasatinib, Gefitinib, and Nilotinib had a higher score in cluster 1 than in cluster 2 (P < 0.05) (Fig. 12). As a result, the above-mentioned 6 drugs were more effective against cluster 2, the hot tumor. These drugs may be used as an adjunctive treatment option to reverse the poor outcome of immunotherapy resistance.

Discussion
For melanoma, the most lethal skin cancer, benefits for patients have vastly improved since the introduction of targeted treatments and immunotherapy [28][29][30] . However, the problems of inadequate response and drug resistance to BRAF/MEK and CTLA-4/PD-1 inhibitors are becoming increasingly tricky 31,32 . Necroptosis, as a form of immunogenic cell death, can overcome apoptosis resistance and effectively facilitate anti-cancer immunity 33,34 . Some lncRNAs have been found to have tumor-suppressor functions, which can eliminate adverse effects of immunotherapy resistance for cancer. On the other hand, some lncRNAs were known to help induce immune system evasion for cancer cells 24,25 . Therefore, the based targets of necroptosis-related lncRNAs would assist with immunotherapy in melanoma as a potential strategy to evade therapy failure.
In the present study, USP30-AS1, LINC01711, LINC00520, NRIR, BASP1-AS1, and LINC02178 were selected as the final components of the risk model. The Sankey diagram shows that only LINC01711 was negatively correlated with BCL2 and AXL among the hub lncRNAs. Silencing of BCL2 in melanoma cells inhibits the polarization of M2 macrophages, promoting anti-cancer immunity and interfering with tumor progression 35 . AXL was prominently expressed in the sentinel lymph nodes of melanoma, suggesting a correlation with metastasis 36 . Overexpression of both genes would affect the outcome of anti-PD-1 immunotherapy 37,38 . LINC02178 has been shown to be associated with lung adenocarcinoma and endometrial cancer, which may serve as a specific prognostic marker for cancer 39,40 . NRIR, firstly found to be associated with melanoma by our study, was positively correlated with CYLD, MLKL, STAT1, and TNFSF10. These genes interfere with melanoma cell progression, epigenetic alterations, the microenvironment of tumor cells, and the effect on immunotherapy 16,41,42 . Besides above 3 lncRNAs, the remaining hub lncRNAs have been linked to the development of melanoma. Existing studies suggest that USP30-AS1 and BASP1-AS1 may be interfere in melanoma cancer cell autophagy 43,44 . Multiple researches have demonstrated that cancer cells with high expression levels of RIPK3 and MLKL are more likely to undergo necroptosis and trigger immune responses that fight tumors 11,45 . USP30-AS1 was highly correlated with RIPK3 and MLKL in the present investigation, however the risk score decreased with increasing USP30-AS1 expression. However, the specific regulatory pathways of USP30-AS1 interfering with RIPK3 and MLKL in melanoma still need further investigation. LINC00520 may interfere with melanoma progression through the miR125b-5p/EIF5A2 axis 46 . LINC00520 expression was positively correlated with AIFM1, and high-risk scores, suggesting a poor outcome. However, the relevance of AIFM1 in melanoma remains uncertain. These obtained hub lncRNAs may provide a foundation for research addressing melanoma treatment failure.
Immune checkpoint inhibitors have transformed the treatment of patients with advanced melanoma, providing many patients with the opportunity to extend their survival 47 . However, many patients do not respond to treatment, and others experience varying degrees of adverse drug effects. Therefore, differentiating between cool and hot tumors by evaluating the immune microenvironment helps better direct immunotherapy interventions 48 . Necroptosis, a potentially novel mechanism of immune cell death, releases large amounts of immunogenic DAMPs. Stimulation of DAMPs may alter the immune resistance of the tumor, transforming it from a cold to a hot tumor. For these reason, necroptosis activation in cancer cells may boost immunotherapy's effectiveness 34 . RNA methylation, a type of RNA modification, happens in both mRNA and lncRNA. Although m7G-related tumor research is still in its infancy, current studies indicate that it interfered with cancer development and can regulate the immune environment of cancer cells by regulating the expression of genes 50,51 . N7-methylguanosine (m7G) has been associated with modifications of a range of genes, particularly METTL1 and WDR4 52 . Our results show that the typical m7G methylation-related genes, METTL1, and WDR4, are more highly expressed in high-risk and cold tumors compared to the other two subgroups. This fantastic result suggests that inhibiting m7G methylation modifications may promote necroptosis in cancer cells and stimulate anti-tumor immunity. However, the concrete mechanisms of these relationships are unclear, which will facilitate future studies of m7G methylation of RNA in melanoma.
With the developing knowledge of the biology underlying melanoma, the combination therapeutic strategy for targeted genetic inhibitors with immune checkpoint inhibitors can significantly benefit advanced patients 29,31 . However, the rising incidence of drug resistance and complications has led to treatment failure and increased patient burden. In the present study, the sensitivity of some common adjuvant melanoma drugs not only differed between high-and low-risk groups but also hot and cold tumors. Therefore, drugs with significant sensitivity may modify the immune microenvironment of the tumor and reduce the adverse effects of immune escape as a complementary option in order to maximize the efficacy of immunotherapy and the survival time for metastatic melanoma patients as much as possible.
Although we attempted to construct the model by increasing the number of random groupings and crossvalidation to make the risk model more accurate, there are still limitations. Above all, as a retrospective study, despite the relatively large sample size and the internal validation, external validation of other databases is still lacking. The reason is related to other datasets with smaller samples, and less complete information is less suitable for validation. Second, there is a dearth of experimental analysis to support the results of bioinformatics analyses. The predicted results for drug sensitivity are based on a pan-cancer public database. Hence, multicenter trial clinical efficacy data and prospective studies are required to demonstrate its viability. Future challenges include a deeper exploration of the mechanisms behind necroptosis-related lncRNAs in melanoma is needed.

Conclusion
Our research successfully predicted the necroptosis-related hub lncRNAs in melanoma, and these lncRNAs provide novel potential targets to eliminate the adverse effects of immunotherapy resistance in melanoma. Furthermore, identifying hot and cold tumors and providing personalized treatment plans that improve the immune microenvironment would provide more significant benefits for advanced patients. Therefore, we consider that necroptosis-related hub lncRNAs have a significant role in aiding earlier detection and improving poor treatment outcomes in malignant melanoma.  www.nature.com/scientificreports/ Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.